#ifndef _VECTORS_H_
#define _VECTORS_H_
#define _CRT_NONSTDC_NO_WARNINGS
//
// originally implemented by Justin Legakis
//

#include <iostream>



#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>

using namespace std;
class Matrix;

// ====================================================================
// ====================================================================

class Vec2f {

public:

	// CONSTRUCTORS & DESTRUCTOR
	Vec2f() { data[0] = data[1] = 0; }
	Vec2f(const Vec2f &V) {
		data[0] = V.data[0];
		data[1] = V.data[1];
	}
	Vec2f(float d0, float d1) {
		data[0] = d0;
		data[1] = d1;
	}
	Vec2f(const Vec2f &V1, const Vec2f &V2) {
		data[0] = V1.data[0] - V2.data[0];
		data[1] = V1.data[1] - V2.data[1];
	}
	~Vec2f() { }

	// ACCESSORS
	void Get(float &d0, float &d1) const {
		d0 = data[0];
		d1 = data[1];
	}
	float operator[](int i) const {
		assert(i >= 0 && i < 2);
		return data[i];
	}
	float x() const { return data[0]; }
	float y() const { return data[1]; }
	float Length() const {
		float l = (float)sqrt(data[0] * data[0] + data[1] * data[1]);
		return l;
	}

	// MODIFIERS
	void Set(float d0, float d1) {
		data[0] = d0;
		data[1] = d1;
	}
	void Scale(float d0, float d1) {
		data[0] *= d0;
		data[1] *= d1;
	}
	void Divide(float d0, float d1) {
		data[0] /= d0;
		data[1] /= d1;
	}
	void Negate() {
		data[0] = -data[0];
		data[1] = -data[1];
	}

	// OVERLOADED OPERATORS
	Vec2f& operator=(const Vec2f &V) {
		data[0] = V.data[0];
		data[1] = V.data[1];
		return *this;
	}
	int operator==(const Vec2f &V) const {
		return ((data[0] == V.data[0]) &&
			(data[1] == V.data[1]));
	}
	int operator!=(const Vec2f &V) {
		return ((data[0] != V.data[0]) ||
			(data[1] != V.data[1]));
	}
	Vec2f& operator+=(const Vec2f &V) {
		data[0] += V.data[0];
		data[1] += V.data[1];
		return *this;
	}
	Vec2f& operator-=(const Vec2f &V) {
		data[0] -= V.data[0];
		data[1] -= V.data[1];
		return *this;
	}
	Vec2f& operator*=(float f) {
		data[0] *= f;
		data[1] *= f;
		return *this;
	}
	Vec2f& operator/=(float f) {
		data[0] /= f;
		data[1] /= f;
		return *this;
	}

	// OPERATIONS
	float Dot2(const Vec2f &V) const {
		return data[0] * V.data[0] + data[1] * V.data[1];
	}

	// STATIC OPERATIONS
	static void Add(Vec2f &a, const Vec2f &b, const Vec2f &c) {
		a.data[0] = b.data[0] + c.data[0];
		a.data[1] = b.data[1] + c.data[1];
	}
	static void Sub(Vec2f &a, const Vec2f &b, const Vec2f &c) {
		a.data[0] = b.data[0] - c.data[0];
		a.data[1] = b.data[1] - c.data[1];
	}
	static void CopyScale(Vec2f &a, const Vec2f &b, float c) {
		a.data[0] = b.data[0] * c;
		a.data[1] = b.data[1] * c;
	}
	static void AddScale(Vec2f &a, const Vec2f &b, const Vec2f &c, float d) {
		a.data[0] = b.data[0] + c.data[0] * d;
		a.data[1] = b.data[1] + c.data[1] * d;
	}
	static void Average(Vec2f &a, const Vec2f &b, const Vec2f &c) {
		a.data[0] = (b.data[0] + c.data[0]) * 0.5f;
		a.data[1] = (b.data[1] + c.data[1]) * 0.5f;
	}
	static void WeightedSum(Vec2f &a, const Vec2f &b, float c, const Vec2f &d, float e) {
		a.data[0] = b.data[0] * c + d.data[0] * e;
		a.data[1] = b.data[1] * c + d.data[1] * e;
	}

	// INPUT / OUTPUT
	void Write(FILE *F = stdout) const {
		fprintf(F, "%f %f\n", data[0], data[1]);
	}

private:

	// REPRESENTATION
	float		data[2];

};

// ====================================================================
// ====================================================================

class Vec3f {

public:

	// CONSTRUCTORS & DESTRUCTOR
	Vec3f() { data[0] = data[1] = data[2] = 0; }
	Vec3f(const Vec3f &V) {
		data[0] = V.data[0];
		data[1] = V.data[1];
		data[2] = V.data[2];
	}
	Vec3f(float d0, float d1, float d2) {
		data[0] = d0;
		data[1] = d1;
		data[2] = d2;
	}
	Vec3f(const Vec3f &V1, const Vec3f &V2) {
		data[0] = V1.data[0] - V2.data[0];
		data[1] = V1.data[1] - V2.data[1];
		data[2] = V1.data[2] - V2.data[2];
	}
	~Vec3f() { }

	// ACCESSORS
	void Get(float &d0, float &d1, float &d2) const {
		d0 = data[0];
		d1 = data[1];
		d2 = data[2];
	}
	float operator[](int i) const {
		assert(i >= 0 && i < 3);
		return data[i];
	}
	float x() const { return data[0]; }
	float y() const { return data[1]; }
	float z() const { return data[2]; }
	float r() const { return data[0]; }
	float g() const { return data[1]; }
	float b() const { return data[2]; }
	float Length() const {
		float l = (float)sqrt(data[0] * data[0] +
			data[1] * data[1] +
			data[2] * data[2]);
		return l;
	}

	// MODIFIERS
	void Set(float d0, float d1, float d2) {
		data[0] = d0;
		data[1] = d1;
		data[2] = d2;
	}
	void Scale(float d0, float d1, float d2) {
		data[0] *= d0;
		data[1] *= d1;
		data[2] *= d2;
	}
	void Divide(float d0, float d1, float d2) {
		data[0] /= d0;
		data[1] /= d1;
		data[2] /= d2;
	}
	void Normalize() {
		float l = Length();
		if (l > 0) {
			data[0] /= l;
			data[1] /= l;
			data[2] /= l;
		}
	}
	void Negate() {
		data[0] = -data[0];
		data[1] = -data[1];
		data[2] = -data[2];
	}
	void Clamp(float low = 0, float high = 1) {
		if (data[0] < low) data[0] = low;  if (data[0] > high) data[0] = high;
		if (data[1] < low) data[1] = low;  if (data[1] > high) data[1] = high;
		if (data[2] < low) data[2] = low;  if (data[2] > high) data[2] = high;
	}

	// OVERLOADED OPERATORS
	Vec3f& operator=(const Vec3f &V) {
		data[0] = V.data[0];
		data[1] = V.data[1];
		data[2] = V.data[2];
		return *this;
	}
	int operator==(const Vec3f &V) {
		return ((data[0] == V.data[0]) &&
			(data[1] == V.data[1]) &&
			(data[2] == V.data[2]));
	}
	int operator!=(const Vec3f &V) {
		return ((data[0] != V.data[0]) ||
			(data[1] != V.data[1]) ||
			(data[2] != V.data[2]));
	}
	Vec3f& operator+=(const Vec3f &V) {
		data[0] += V.data[0];
		data[1] += V.data[1];
		data[2] += V.data[2];
		return *this;
	}
	Vec3f& operator-=(const Vec3f &V) {
		data[0] -= V.data[0];
		data[1] -= V.data[1];
		data[2] -= V.data[2];
		return *this;
	}
	Vec3f& operator*=(int i) {
		data[0] = float(data[0] * i);
		data[1] = float(data[1] * i);
		data[2] = float(data[2] * i);
		return *this;
	}
	Vec3f& operator*=(float f) {
		data[0] *= f;
		data[1] *= f;
		data[2] *= f;
		return *this;
	}
	Vec3f& operator/=(int i) {
		data[0] = float(data[0] / i);
		data[1] = float(data[1] / i);
		data[2] = float(data[2] / i);
		return *this;
	}
	Vec3f& operator/=(float f) {
		data[0] /= f;
		data[1] /= f;
		data[2] /= f;
		return *this;
	}


	friend Vec3f operator+(const Vec3f &v1, const Vec3f &v2) {
		Vec3f v3; Add(v3, v1, v2); return v3;
	}
	friend Vec3f operator-(const Vec3f &v1, const Vec3f &v2) {
		Vec3f v3; Sub(v3, v1, v2); return v3;
	}
	friend Vec3f operator*(const Vec3f &v1, float f) {
		Vec3f v2; CopyScale(v2, v1, f); return v2;
	}
	friend Vec3f operator*(float f, const Vec3f &v1) {
		Vec3f v2; CopyScale(v2, v1, f); return v2;
	}
	friend Vec3f operator*(const Vec3f &v1, const Vec3f &v2) {
		Vec3f v3; Mult(v3, v1, v2); return v3;
	}


	// OPERATIONS
	float Dot3(const Vec3f &V) const {
		return data[0] * V.data[0] +
			data[1] * V.data[1] +
			data[2] * V.data[2];
	}

	// STATIC OPERATIONS
	static void Add(Vec3f &a, const Vec3f &b, const Vec3f &c) {
		a.data[0] = b.data[0] + c.data[0];
		a.data[1] = b.data[1] + c.data[1];
		a.data[2] = b.data[2] + c.data[2];
	}
	static void Sub(Vec3f &a, const Vec3f &b, const Vec3f &c) {
		a.data[0] = b.data[0] - c.data[0];
		a.data[1] = b.data[1] - c.data[1];
		a.data[2] = b.data[2] - c.data[2];
	}
	static void Mult(Vec3f &a, const Vec3f &b, const Vec3f &c) {
		a.data[0] = b.data[0] * c.data[0];
		a.data[1] = b.data[1] * c.data[1];
		a.data[2] = b.data[2] * c.data[2];
	}
	static void CopyScale(Vec3f &a, const Vec3f &b, float c) {
		a.data[0] = b.data[0] * c;
		a.data[1] = b.data[1] * c;
		a.data[2] = b.data[2] * c;
	}
	static void AddScale(Vec3f &a, const Vec3f &b, const Vec3f &c, float d) {
		a.data[0] = b.data[0] + c.data[0] * d;
		a.data[1] = b.data[1] + c.data[1] * d;
		a.data[2] = b.data[2] + c.data[2] * d;
	}
	static void Average(Vec3f &a, const Vec3f &b, const Vec3f &c) {
		a.data[0] = (b.data[0] + c.data[0]) * 0.5f;
		a.data[1] = (b.data[1] + c.data[1]) * 0.5f;
		a.data[2] = (b.data[2] + c.data[2]) * 0.5f;
	}
	static void WeightedSum(Vec3f &a, const Vec3f &b, float c, const Vec3f &d, float e) {
		a.data[0] = b.data[0] * c + d.data[0] * e;
		a.data[1] = b.data[1] * c + d.data[1] * e;
		a.data[2] = b.data[2] * c + d.data[2] * e;
	}
	static void Cross3(Vec3f &c, const Vec3f &v1, const Vec3f &v2) {
		float x = v1.data[1] * v2.data[2] - v1.data[2] * v2.data[1];
		float y = v1.data[2] * v2.data[0] - v1.data[0] * v2.data[2];
		float z = v1.data[0] * v2.data[1] - v1.data[1] * v2.data[0];
		c.data[0] = x; c.data[1] = y; c.data[2] = z;
	}

	static void Min(Vec3f &a, const Vec3f &b, const Vec3f &c) {
		a.data[0] = (b.data[0] < c.data[0]) ? b.data[0] : c.data[0];
		a.data[1] = (b.data[1] < c.data[1]) ? b.data[1] : c.data[1];
		a.data[2] = (b.data[2] < c.data[2]) ? b.data[2] : c.data[2];
	}
	static void Max(Vec3f &a, const Vec3f &b, const Vec3f &c) {
		a.data[0] = (b.data[0] > c.data[0]) ? b.data[0] : c.data[0];
		a.data[1] = (b.data[1] > c.data[1]) ? b.data[1] : c.data[1];
		a.data[2] = (b.data[2] > c.data[2]) ? b.data[2] : c.data[2];
	}

	// INPUT / OUTPUT
	void Write(FILE *F = stdout) const {
		fprintf(F, "%f %f %f\n", data[0], data[1], data[2]);
	}

private:

	friend class Matrix;

	// REPRESENTATION
	float		data[3];

};

// ====================================================================
// ====================================================================

class Vec4f {

public:

	// CONSTRUCTORS & DESTRUCTOR
	Vec4f() { data[0] = data[1] = data[2] = data[3] = 0; }
	Vec4f(const Vec4f &V) {
		data[0] = V.data[0];
		data[1] = V.data[1];
		data[2] = V.data[2];
		data[3] = V.data[3];
	}
	Vec4f(float d0, float d1, float d2, float d3) {
		data[0] = d0;
		data[1] = d1;
		data[2] = d2;
		data[3] = d3;
	}
	Vec4f(const Vec3f &V, float w) {
		data[0] = V.x();
		data[1] = V.y();
		data[2] = V.z();
		data[3] = w;
	}
	Vec4f(const Vec4f &V1, const Vec4f &V2) {
		data[0] = V1.data[0] - V2.data[0];
		data[1] = V1.data[1] - V2.data[1];
		data[2] = V1.data[2] - V2.data[2];
		data[3] = V1.data[3] - V2.data[3];
	}
	~Vec4f() { }

	// ACCESSORS
	void Get(float &d0, float &d1, float &d2, float &d3) const {
		d0 = data[0];
		d1 = data[1];
		d2 = data[2];
		d3 = data[3];
	}
	float operator[](int i) const {
		assert(i >= 0 && i < 4);
		return data[i];
	}
	float x() const { return data[0]; }
	float y() const { return data[1]; }
	float z() const { return data[2]; }
	float w() const { return data[3]; }
	float r() const { return data[0]; }
	float g() const { return data[1]; }
	float b() const { return data[2]; }
	float a() const { return data[3]; }
	float Length() const {
		float l = (float)sqrt(data[0] * data[0] +
			data[1] * data[1] +
			data[2] * data[2] +
			data[3] * data[3]);
		return l;
	}

	// MODIFIERS
	void Set(float d0, float d1, float d2, float d3) {
		data[0] = d0;
		data[1] = d1;
		data[2] = d2;
		data[3] = d3;
	}
	void Scale(float d0, float d1, float d2, float d3) {
		data[0] *= d0;
		data[1] *= d1;
		data[2] *= d2;
		data[3] *= d3;
	}
	void Divide(float d0, float d1, float d2, float d3) {
		data[0] /= d0;
		data[1] /= d1;
		data[2] /= d2;
		data[3] /= d3;
	}
	void Negate() {
		data[0] = -data[0];
		data[1] = -data[1];
		data[2] = -data[2];
		data[3] = -data[3];
	}
	void Normalize() {
		float l = Length();
		if (l > 0) {
			data[0] /= l;
			data[1] /= l;
			data[2] /= l;
		}
	}
	void DivideByW() {
		if (data[3] != 0) {
			data[0] /= data[3];
			data[1] /= data[3];
			data[2] /= data[3];
		}
		else {
			data[0] = data[1] = data[2] = 0;
		}
		data[3] = 1;
	}

	// OVERLOADED OPERATORS
	Vec4f& operator=(const Vec4f &V) {
		data[0] = V.data[0];
		data[1] = V.data[1];
		data[2] = V.data[2];
		data[3] = V.data[3];
		return *this;
	}
	int operator==(const Vec4f &V) const {
		return ((data[0] == V.data[0]) &&
			(data[1] == V.data[1]) &&
			(data[2] == V.data[2]) &&
			(data[3] == V.data[3]));
	}
	int operator!=(const Vec4f &V) const {
		return ((data[0] != V.data[0]) ||
			(data[1] != V.data[1]) ||
			(data[2] != V.data[2]) ||
			(data[3] != V.data[3]));
	}
	Vec4f& operator+=(const Vec4f &V) {
		data[0] += V.data[0];
		data[1] += V.data[1];
		data[2] += V.data[2];
		data[3] += V.data[3];
		return *this;
	}
	Vec4f& operator-=(const Vec4f &V) {
		data[0] -= V.data[0];
		data[1] -= V.data[1];
		data[2] -= V.data[2];
		data[3] -= V.data[3];
		return *this;
	}
	Vec4f& operator*=(float f) {
		data[0] *= f;
		data[1] *= f;
		data[2] *= f;
		data[3] *= f;
		return *this;
	}
	Vec4f& operator/=(float f) {
		data[0] /= f;
		data[1] /= f;
		data[2] /= f;
		data[3] /= f;
		return *this;
	}

	// OPERATIONS
	float Dot2(const Vec4f &V) const {
		return data[0] * V.data[0] +
			data[1] * V.data[1];
	}
	float Dot3(const Vec4f &V) const {
		return data[0] * V.data[0] +
			data[1] * V.data[1] +
			data[2] * V.data[2];
	}
	float Dot4(const Vec4f &V) const {
		return data[0] * V.data[0] +
			data[1] * V.data[1] +
			data[2] * V.data[2] +
			data[3] * V.data[3];
	}

	// STATIC OPERATIONS
	static void Add(Vec4f &a, const Vec4f &b, const Vec4f &c) {
		a.data[0] = b.data[0] + c.data[0];
		a.data[1] = b.data[1] + c.data[1];
		a.data[2] = b.data[2] + c.data[2];
		a.data[3] = b.data[3] + c.data[3];
	}
	static void Sub(Vec4f &a, const Vec4f &b, const Vec4f &c) {
		a.data[0] = b.data[0] - c.data[0];
		a.data[1] = b.data[1] - c.data[1];
		a.data[2] = b.data[2] - c.data[2];
		a.data[3] = b.data[3] - c.data[3];
	}
	static void CopyScale(Vec4f &a, const Vec4f &b, float c) {
		a.data[0] = b.data[0] * c;
		a.data[1] = b.data[1] * c;
		a.data[2] = b.data[2] * c;
		a.data[3] = b.data[3] * c;
	}
	static void AddScale(Vec4f &a, const Vec4f &b, const Vec4f &c, float d) {
		a.data[0] = b.data[0] + c.data[0] * d;
		a.data[1] = b.data[1] + c.data[1] * d;
		a.data[2] = b.data[2] + c.data[2] * d;
		a.data[3] = b.data[3] + c.data[3] * d;
	}
	static void Average(Vec4f &a, const Vec4f &b, const Vec4f &c) {
		a.data[0] = (b.data[0] + c.data[0]) * 0.5f;
		a.data[1] = (b.data[1] + c.data[1]) * 0.5f;
		a.data[2] = (b.data[2] + c.data[2]) * 0.5f;
		a.data[3] = (b.data[3] + c.data[3]) * 0.5f;
	}
	static void WeightedSum(Vec4f &a, const Vec4f &b, float c, const Vec4f &d, float e) {
		a.data[0] = b.data[0] * c + d.data[0] * e;
		a.data[1] = b.data[1] * c + d.data[1] * e;
		a.data[2] = b.data[2] * c + d.data[2] * e;
		a.data[3] = b.data[3] * c + d.data[3] * e;
	}
	static void Cross3(Vec4f &c, const Vec4f &v1, const Vec4f &v2) {
		float x = v1.data[1] * v2.data[2] - v1.data[2] * v2.data[1];
		float y = v1.data[2] * v2.data[0] - v1.data[0] * v2.data[2];
		float z = v1.data[0] * v2.data[1] - v1.data[1] * v2.data[0];
		c.data[0] = x; c.data[1] = y; c.data[2] = z;
	}

	// INPUT / OUTPUT
	void Write(FILE *F = stdout) const {
		fprintf(F, "%f %f %f %f\n", data[0], data[1], data[2], data[3]);
	}

private:

	friend class Matrix;

	// REPRESENTATION
	float		data[4];

};

inline ostream &operator<<(ostream &os, const Vec3f &v) {
	os << "Vec3f <" << v.x() << ", " << v.y() << ", " << v.z() << ">";
	return os;
}

// ====================================================================
// ====================================================================

#endif

